{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0d3f0e73",
   "metadata": {},
   "source": [
    "# 机器读心术之文本挖掘与自然语言处理第5课书面作业\n",
    "学号：207402  \n",
    "\n",
    "**书面作业**  \n",
    "1. 在上周作业的基础上，体验HMM第二个基本问题的维特比算法，计算上周幻灯片第34页“掷骰子HMM”中可观测序列为“1635273524”，反求出最大可能的隐藏状态序列，可使用任何编程语言或手工计算。  \n",
    "2. 假设语料库为  \n",
    "研究生物很有意思  \n",
    "他大学时代是研究生物的  \n",
    "生物专业是他的首选目标  \n",
    "他是研究生  \n",
    "\n",
    "（1）试直接用肉眼观察加注分词标注（BMES），然后基于此语料建立相应的HMM，求解出全套模型参数  \n",
    "（2）【可选】如果只有原始语料，没有分词标注，可以建立HMM吗？如果能建立，训练过程是怎样的？  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba9f5665",
   "metadata": {},
   "source": [
    "## 作业1\n",
    "在上周作业的基础上，体验HMM第二个基本问题的维特比算法，计算上周幻灯片第34页“掷骰子HMM”中可观测序列为“1635273524”，反求出最大可能的隐藏状态序列，可使用任何编程语言或手工计算。  \n",
    "**解：**  \n",
    "根据宗庆华书中维特比算法定义如下："
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e37b167",
   "metadata": {},
   "source": [
    "维特比变量$\\delta_t(x)$是在时间$t$时，HMM沿着某一条路径到达状态$s_i$，并输出观察序列$O_1O_2...O_t$的最大概率：\n",
    "$$\n",
    "\\delta_t(i)=\\underset{q_1,q_2,\\ldots,q_{t-1}}{\\operatorname{max}}P(q_1,q_2,\\ldots,q_t=s_i,O_1O_2\\ldots O_t|\\mu) \\tag{6.21}\n",
    "$$\n",
    "同时设置另外一个变量$\\psi_t(i)$用于路径记忆，让$\\psi_t(i)$记录该路径上状态$s_i$的前一个（在时间$t-1$）状态。  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f859fac",
   "metadata": {},
   "source": [
    "**步1** 初始化：  \n",
    "$$\n",
    "\\begin{align*}\n",
    "\\delta_1(i)&=\\pi_i b_i(O_1) \\quad,\\quad 1 \\le i \\le N  \\\\\n",
    "\\psi_1&=0\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cebf09af",
   "metadata": {},
   "source": [
    "**步2** 归纳计算：\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\delta_t(j)&=\\underset{1\\le i\\le N}{\\operatorname{max}}[\\delta_{t-1}(i)\\cdot a_{ij}]\\cdot b_j(O_t)\\quad ,\\quad 2\\le t \\le T;1\\le j\\le N \\\\\n",
    "\\psi_t(j)&=\\underset{1\\le i\\le N}{\\operatorname{arg \\, max}}[\\delta_{t-1}(i)\\cdot a_{ij}]\\cdot b_j(O_t)\\quad ,\\quad 2\\le t\\le T;1\\le i\\le N\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edb4d7fa",
   "metadata": {},
   "source": [
    "上面关于$\\psi_t(j)$的公式应该不对，后面不需要乘以$b_j(O_t)$:\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\psi_t(j)&=\\underset{1\\le i\\le N}{\\operatorname{arg \\, max}}[\\delta_{t-1}(i)\\cdot a_{ij}]\\quad ,\\quad 2\\le t\\le T;1\\le i\\le N\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e6d773c",
   "metadata": {},
   "source": [
    "**步3** 终结：\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\hat{Q}_T&=\\underset{1\\le i\\le N}{\\operatorname{arg\\,max}}[\\delta_T(i)] \\\\\n",
    "\\hat{P}(\\hat{Q}_T)&=\\underset{1\\le i\\le N}{\\operatorname{max}}[\\delta_T(i)]\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a588ad0",
   "metadata": {},
   "source": [
    "**步4** 路径（状态序列）回溯：\n",
    "$$\n",
    "\\hat{q}_t=\\psi_{t+1}(\\hat{q}_{t+1})\\quad ,\\quad t=T-1,T-2,\\ldots,1\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "95090617",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# 状态转换概率矩阵定义如下：\n",
    "a=np.array(\n",
    "    [\n",
    "        [1/3.,1/3.,1/3.],\n",
    "        [1/3.,1/3.,1/3.],\n",
    "        [1/3.,1/3.,1/3.]\n",
    "    ]\n",
    ")\n",
    "\n",
    "# 发射概率矩阵定义如下：\n",
    "b=np.array(\n",
    "    [\n",
    "        [1/6,1/6,1/6,1/6,1/6,1/6,  0,  0],\n",
    "        [1/4,1/4,1/4,1/4,  0,  0,  0,  0],\n",
    "        [1/8,1/8,1/8,1/8,1/8,1/8,1/8,1/8]\n",
    "    ]\n",
    ")\n",
    "\n",
    "#初始概率分布定义如下：\n",
    "pi=np.array([1/3,1/3,1/3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "6ae1dc72",
   "metadata": {},
   "outputs": [],
   "source": [
    "def viterbi_proc(a,b,pi,seq,N):\n",
    "    '''\n",
    "    维特比算法\n",
    "    a: 状态转换概率矩阵\n",
    "    b: 发射概率矩阵\n",
    "    pi: 初始概率分布\n",
    "    seq: 观察序列\n",
    "    N: 状态数\n",
    "    '''\n",
    "    T=len(seq)\n",
    "    delta=[]\n",
    "    delta.append(pi*b[:,0])#第1步，初始化\n",
    "    psi=[]\n",
    "    \n",
    "    for t in range(1,T,1):\n",
    "        delta_t=np.zeros(N)\n",
    "        psi_t=np.zeros(N)\n",
    "        for j in range(N):\n",
    "            delta_t_sub_1=delta[-1]\n",
    "            delta_t[j]=np.max(delta_t_sub_1 * a[:,j])*b[j,int(seq[t])]\n",
    "            psi_t[j]=np.argmax(delta_t_sub_1 * a[:,j])\n",
    "        delta.append(delta_t)\n",
    "        psi.append(psi_t)\n",
    "    rpsi=[]\n",
    "    rpsi.append(np.argmax(delta[-1]))\n",
    "    psi.pop()\n",
    "    L=len(psi)\n",
    "    for i in range(L):\n",
    "        t=L-i-1\n",
    "        rpsi.append(int(psi[t][int(rpsi[-1])]))\n",
    "    rpsi.reverse()\n",
    "    return np.max(delta[-1]),rpsi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "14ae0fc8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1.1963376510358935e-12, [1, 2, 1, 0, 1, 2, 1, 0, 0])\n"
     ]
    }
   ],
   "source": [
    "print(viterbi_proc(a,b,pi,'1635273524',3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12aef861",
   "metadata": {},
   "source": [
    "## 作业2\n",
    "\n",
    "假设语料库为：  \n",
    "* 研究生物很有意思  \n",
    "* 他大学时代是研究生物的\n",
    "* 生物专业是他的首选目标\n",
    "* 他是研究生  \n",
    "\n",
    "（1）试直接用肉眼观察加注分词标注（BMES），然后基于此语料建立相应的HMM，求解出全套模型参数  \n",
    "（2）【可选】如果只有原始语料，没有分词标注，可以建立HMM吗？如果能建立，训练过程是怎样的？  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b144d11",
   "metadata": {},
   "source": [
    "**解：**  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2eb3f579",
   "metadata": {},
   "source": [
    "(1)  \n",
    "先进行肉眼标注：  \n",
    "* $\\overset{B}{研}\\overset{E}{究}\\overset{B}{生}\\overset{E}{物}\\overset{S}{很}\\overset{S}{有}\\overset{B}{意}\\overset{E}{思}$  \n",
    "* $\\overset{S}{他}\\overset{B}{大}\\overset{E}{学}\\overset{B}{时}\\overset{E}{代}\\overset{S}{是}\\overset{B}{研}\\overset{E}{究}\\overset{B}{生}\\overset{E}{物}\\overset{S}{的}$  \n",
    "* $\\overset{B}{生}\\overset{E}{物}\\overset{B}{专}\\overset{E}{业}\\overset{S}{是}\\overset{S}{他}\\overset{S}{的}\\overset{B}{首}\\overset{E}{选}\\overset{B}{目}\\overset{E}{标}$  \n",
    "* $\\overset{S}{他}\\overset{S}{是}\\overset{B}{研}\\overset{M}{究}\\overset{E}{生}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc1cf0d0",
   "metadata": {},
   "source": [
    "在本例中，隐藏序列是BMES，共有4个状态即B、M、E、S。观察序列就是上面的中文四句话。我们要求解全套模型参数，即：$\\pi(i),a_{ij},b_j(k)$，其中$1\\le i,j \\le 4$，这里令$s_1$=B, $s_2$=M,$s_3$=E,$s_4$=S。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ca4b711",
   "metadata": {},
   "source": [
    "这里共有汉字21个，因此$1\\le k \\le 21$。这里令：  \n",
    "$$\n",
    "q_1=研, q_2=究, q_3=生, q_4=物, q_5=很, q_6=有, q_7=意 \\\\\n",
    "q_8=思, q_9=他, q_{10}=大, q_{11}=学, q_{12}=时, q_{13}=代, q_{14}=是 \\\\\n",
    "q_{15}=的, q_{16}=专, q_{17}=业, q_{18}=首, q_{19}=选, q_{20}=目, q_{21}=标 \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bec543e0",
   "metadata": {},
   "source": [
    "根据语料库及肉眼标注情况可得：  \n",
    "$$\n",
    "\\pi=\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "\\frac{1}{2} & 0 & 0 & \\frac{1}{2}\\\\\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42569e5c",
   "metadata": {},
   "source": [
    "$$\n",
    "a=\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "0 & \\frac{1}{12} & \\frac{11}{12} & 0\\\\\n",
    "0 & 0 & 1 & 0\\\\\n",
    "\\frac{5}{12} & 0 & 0 & \\frac{4}{12}\\\\\n",
    "\\frac{1}{2} & 0 & 0 & \\frac{2}{5}\\\\\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52104eda",
   "metadata": {},
   "source": [
    "$$\n",
    "b=\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "\\frac{3}{12} & 0 & \\frac{3}{12} & 0 & 0 & 0 & \\frac{1}{12} & 0 & 0 & \\frac{1}{12} & 0 & \\frac{1}{12} & 0 & 0 & 0 & \\frac{1}{12} & 0 & \\frac{1}{12} & 0 & \\frac{1}{12} & 0  \\\\\n",
    "0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\\\\n",
    "0 & \\frac{2}{12} & \\frac{1}{12} & \\frac{3}{12} & 0 & 0 & 0 & \\frac{1}{12} & 0 & 0 & \\frac{1}{12} & 0 & \\frac{1}{12} & 0 & 0 & 0 & \\frac{1}{12} & 0 & \\frac{1}{12} & 0 & \\frac{1}{12}\\\\\n",
    "0 & 0 & 0 & 0 & \\frac{1}{12} & \\frac{1}{12} & 0 & 0 & \\frac{3}{12} & 0 & 0 & 0 & 0 & \\frac{3}{12} & \\frac{2}{12} & 0 & 0 & 0 & 0 & 0 & 0\\\\\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2b2bbf2",
   "metadata": {},
   "source": [
    "(2)  \n",
    "如果只有原始语料，没有分词标注，也可以建立HMM，可以根据初始给定，利用EM算法计算HMM模型参数。根据本例计算如下："
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3a13de6",
   "metadata": {},
   "source": [
    "根据宗庆华书中的向前向后算法：  \n",
    "**步1** 初始化：随机地给参数$\\pi(i),a_{ij},b_j(k)$赋值，使其满足如下条件： "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bdd477e",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\sum_{i=1}^N \\pi(i)&=1 \\\\\n",
    "\\sum_{j=1}^N a_{ij}&=1 \\quad,\\: 1 \\le i \\le N \\\\\n",
    "\\sum_{k=1}^M b_j(k)&=1\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69eebbc1",
   "metadata": {},
   "source": [
    "**步2** EM计算：  \n",
    "**E-步骤** 由模型$\\mu_i$根据公式(6.24)和(6.25)计算期望值$\\epsilon_t(i,j),\\gamma_t(i)$。  \n",
    "**M-步骤** 用E-步骤得到的$\\epsilon_t(i,j),\\gamma_t(i)$，根据公式(6.26)(6.27)(6.28)重新估计参数$\\pi(i),a_{ij},b_j(k)$的值，从而得到模型$\\mu_{i+1}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04649c44",
   "metadata": {},
   "source": [
    "**步3** 循环计算：  \n",
    "令i=i+1，重复执行EM计算，直到$\\pi(i),a_{ij},b_j(k)$收敛。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b63f7230",
   "metadata": {},
   "source": [
    "$$\n",
    "\\alpha_t(i)=P(O_1O_2\\cdots O_t|q_t=s_i,\\mu) \\tag{6.12}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7523a76",
   "metadata": {},
   "source": [
    "$$\n",
    "\\beta_t(i)=P(O_{t+1}O_{t+2}\\cdots O_T|q_t=s_i,\\mu) \\tag{6.15}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4cb16a25",
   "metadata": {},
   "source": [
    "给定HMM $\\mu$和观察序列$O=O_1O_2\\cdots O_T$，在时间$t$位于状态$s_i$，在时间$t+1$位于状态$s_j$的概率$\\epsilon_t(i,j)$为\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\epsilon_t(i,j)&=P(q_t=s_i,q_{t+1}=s_j|O,\\mu)   \\\\\n",
    "&=\\frac{P(q_t=s_i,q_{t+1}=s_j, O|\\mu)}{P(O|\\mu)} \\\\\n",
    "&=\\frac{\\alpha_t(i)a_{ij}b_j(O_{t+1})\\beta_{t+1}(j)}{\\sum_{i=1}^N\\sum_{j=1}^N \\alpha_t(i)a_{ij}b_j(O_{t+1})\\beta_{t+1}(j)} \\tag{6.24}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc45bc89",
   "metadata": {},
   "source": [
    "给定HMM $\\mu$和观察序列$O=O_1O_2\\cdots O_T$，在时间$t$位于状态$s_i$的概率$\\gamma_t(i)$为\n",
    "$$\n",
    "\\gamma_t(i)=\\sum_{j=1}^N \\epsilon_t(i,j) \\tag{6.25}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "507918a5",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\bar{\\pi}_i&=P(q1=si|O,\\mu)=\\gamma_1(i) \\tag{6.26}\\\\\n",
    "\\bar{a}_{ij}&=\\frac{q_i转移到q_j的次数}{q_i转移到任一态的次数 }  \\\\\n",
    "&=\\frac{\\sum_{t=1}^{T-1}\\epsilon_t(i,j)}{\\sum_{t=1}^{T-1}\\gamma_t(i)} \\tag{6.27}\\\\\n",
    "\\bar{b}_j(k)&=\\frac{Q中从状态q_j发出v_k的次数}{Q中到达q_j的次数 } \\\\\n",
    "&=\\frac{\\sum_{t=1}^{T}\\gamma_t(j)\\delta(O_t,v_k)}{\\sum_{t=1}^{T}\\gamma_t(j)} \\tag{6.28}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d72f1d07",
   "metadata": {},
   "source": [
    "实现源代码如下：  \n",
    "先给参数们一个初始值,这里初始值是按均匀分布随意给的。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "be913004",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# 状态转换概率矩阵定义如下：\n",
    "a=np.array(\n",
    "    [\n",
    "        [1/4.,1/4.,1/4.,1/4.],\n",
    "        [1/4.,1/4.,1/4.,1/4.],\n",
    "        [1/4.,1/4.,1/4.,1/4.],\n",
    "        [1/4.,1/4.,1/4.,1/4.]\n",
    "    ]\n",
    ")\n",
    "\n",
    "# 发射概率矩阵定义如下：\n",
    "b=np.array(\n",
    "    [\n",
    "        [1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.],\n",
    "        [1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.],\n",
    "        [1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.],\n",
    "        [1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.]\n",
    "    ]\n",
    ")\n",
    "\n",
    "#初始概率分布定义如下：\n",
    "pi=np.array([1/4.,1/4.,1/4.,1/4.])\n",
    "\n",
    "#定义观察序列\n",
    "q={'研': 0, '究': 1, '生': 2, '物': 3,'很': 4,'有': 5,'意': 6,\n",
    "   '思': 7, '他': 8, '大': 9, '学':10,'时':11,'代':12,'是':13,\n",
    "   '的':14, '专':15, '业':16, '首':17,'选':18,'目':19,'标':20}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "e11746d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "class HMM:\n",
    "    def __init__(self,A,B,pi):\n",
    "        self.A =A\n",
    "        self.B =B\n",
    "        self.pi =pi\n",
    "    def _forward(self,obs_seq):\n",
    "        # 取A = N x N\n",
    "        N = self.A.shape[0]\n",
    "        T = len(obs_seq)\n",
    "\n",
    "        F = np.zeros((N,T))\n",
    "\n",
    "        # alpha = pi*b\n",
    "        F[:,0] = self.pi *self.B[:,obs_seq[0]]\n",
    "\n",
    "        for t in range(1,T):\n",
    "            for n in range(N):\n",
    "                # 计算第t时，第n个状态的前向概率\n",
    "                F[n,t] = np.dot(F[:,t-1], (self.A[:,n])) * self.B[n, obs_seq[t]]\n",
    "        return F\n",
    "    def _backward(self,obs_seq):\n",
    "        N = self.A.shape[0]\n",
    "        T = len(obs_seq)\n",
    "\n",
    "        X = np.zeros((N,T))\n",
    "        # 表示X矩阵的最后一列\n",
    "        X[:,-1:] =1\n",
    "\n",
    "        for t in reversed(range(T-1)):\n",
    "            for n in range(N):\n",
    "                # 边权值为a_ji\n",
    "                X[n,t] = np.sum(X[:,t+1] * self.A[n,:] * self.B[:,obs_seq[t+1]])\n",
    "\n",
    "        return X\n",
    "\n",
    "    def forward_backward(self, observations, criterion=0.05):\n",
    "        n_states = self.A.shape[0]\n",
    "        # 观察序列的长度T\n",
    "        n_samples = len(observations)\n",
    "\n",
    "        done = False\n",
    "        while not done:\n",
    "            # alpha_t(i) = P(o_1,o_2,...,o_t,q_t = s_i | hmm)\n",
    "            # Initialize alpha\n",
    "            # 获得所有前向传播节点值 alpha_t(i)\n",
    "            alpha = self._forward(observations)\n",
    "\n",
    "            # beta_t(i) = P(o_t+1,o_t+2,...,o_T | q_t = s_i , hmm)\n",
    "            # Initialize beta\n",
    "            # 获得所有后向传播节点值 beta_t(i)\n",
    "            beta = self._backward(observations)\n",
    "\n",
    "            # 计算 xi_t(i,j) -> xi(i,j,t)\n",
    "            xi = np.zeros((n_states, n_states, n_samples - 1))\n",
    "            # 在每个时刻\n",
    "            for t in range(n_samples - 1):\n",
    "                # 计算P(O | hmm)\n",
    "                denom = np.sum(alpha[:, -1])\n",
    "                for i in range(n_states):\n",
    "                    # numer[1,:] = 行向量，alpha[i,t]=实数，slef.A[i,:] = 行向量\n",
    "                    # self.B[:,observations[t+1]].T = 行向量,beta[:,t+1].T = 行向量\n",
    "                    numer = alpha[i, t] * self.A[i, :] * self.B[:, observations[t + 1]].T * beta[:, t + 1].T\n",
    "                    xi[i, :, t] = numer / denom\n",
    "\n",
    "                # 计算gamma_t(i) 就是对j进行求和\n",
    "                gamma = np.sum(xi, axis=1)\n",
    "                # need final gamma elements for new B\n",
    "                prod = (alpha[:, n_samples - 1] * beta[:, n_samples - 1]).reshape((-1, 1))\n",
    "                # 合并T时刻的节点\n",
    "                gamma = np.hstack((gamma, prod / sum(prod)))\n",
    "                # 列向量\n",
    "                newpi = gamma[:, 0]\n",
    "                newA = np.sum(xi, 2) / np.sum(gamma[:, :-1], axis=1).reshape((-1, 1))\n",
    "                newB = np.copy(self.B)\n",
    "\n",
    "                # 观测状态数\n",
    "                num_levels = self.B.shape[1]\n",
    "                sumgamma = np.sum(gamma, axis=1)\n",
    "                for lev in range(num_levels):\n",
    "                    mask = observations == lev\n",
    "                    newB[:, lev] = np.sum(gamma[:, mask], axis=1) / sumgamma\n",
    "\n",
    "                if np.max(np.abs(self.pi - newpi)) < criterion and \\\n",
    "                                np.max(np.abs(self.A - newA)) < criterion and \\\n",
    "                                np.max(np.abs(self.B - newB)) < criterion:\n",
    "                    done = True\n",
    "                self.A[:], self.B[:], self.pi[:] = newA, newB, newpi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "0cc83216",
   "metadata": {},
   "outputs": [],
   "source": [
    "def observation_list(ob):\n",
    "    ol=[]\n",
    "    T=len(ob)\n",
    "    for i in range(T):\n",
    "        ol.append(q[ob[i]])\n",
    "    return np.array(ol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "ccf143db",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A: \n",
      " [[0.25 0.25 0.25 0.25]\n",
      " [0.25 0.25 0.25 0.25]\n",
      " [0.25 0.25 0.25 0.25]\n",
      " [0.25 0.25 0.25 0.25]]\n",
      "B: \n",
      " [[0.08 0.   0.   0.   0.   0.   0.84 0.08 0.   0.   0.   0.   0.   0.\n",
      "  0.   0.   0.   0.   0.   0.   0.  ]\n",
      " [0.08 0.   0.   0.   0.   0.   0.84 0.08 0.   0.   0.   0.   0.   0.\n",
      "  0.   0.   0.   0.   0.   0.   0.  ]\n",
      " [0.08 0.   0.   0.   0.   0.   0.84 0.08 0.   0.   0.   0.   0.   0.\n",
      "  0.   0.   0.   0.   0.   0.   0.  ]\n",
      " [0.08 0.   0.   0.   0.   0.   0.84 0.08 0.   0.   0.   0.   0.   0.\n",
      "  0.   0.   0.   0.   0.   0.   0.  ]]\n",
      "PI: \n",
      " [0.25 0.25 0.25 0.25]\n"
     ]
    }
   ],
   "source": [
    "hmm=HMM(a,b,pi)\n",
    "hmm.forward_backward(observation_list('研究生物很有意思'),criterion=0.01)\n",
    "print('A: \\n',hmm.A)\n",
    "print('B: \\n',hmm.B)\n",
    "print('PI: \\n',hmm.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7c748d0",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
